Atomic-scale insights on the formation of ordered arrays of edge dislocations in Ge/Si(001) films via molecular dynamics simulations

Heteroepitaxial films of Ge on Si(001) are receiving wide attention due to several possible applications in micro- and opto-electronics. Understanding the dynamic behavior of linear defects, such as dislocations, is key. They are unavoidably present in such systems due to the lattice mismatch between the two materials, and can directly influence devices performances. It has been experimentally demonstrated more than fifteen years ago that a suitable choice of the growth parameters allows for the formation of a nicely ordered net of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90^{\circ }$$\end{document}90∘ dislocations at the Ge/Si interface, improving the overall film quality and strain relaxation uniformity. Atomic-scale details on the set of mechanisms leading to such an outcome are however still missing. Here we present a set of classical molecular dynamics simulations shedding light on the full set of microscopic processes driving to the experimentally observed array of linear defects. This includes simple gliding of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$60^{\circ }$$\end{document}60∘ dislocations and vacancy-promoted climbing and gliding. The importance of the particular experimental conditions, involving a low-temperature stage followed by an increase in temperature, is highlighted.

www.nature.com/scientificreports/ demonstrated that thermal treatment in the form of annealing leads to a formation of an ordered network of 90 • (edge) dislocation segments at the interface, improving film quality and strain relaxation uniformity 28,29 . Edge misfit dislocations (often called Lomer dislocations) are characterized by Burgers vectors � b = ±a/2[110] or � b = ±a/2 [110] lying in the interface plane. Two different questions are still under discussion: (i) How are 90 • dislocations formed? (ii) How is lateral ordering of dislocation at the interface reached?
A partial answer to the first question was supplied as different possible mechanisms leading to the formation of edge dislocations at the interface were proposed 30,31 . Experiments and theoretical calculations [31][32][33][34] have highlighted the importance of one of them, i.e. preferential nucleation of a "complementary" 35 60 • dislocation: the stress field induced in the proximity of the free surface by the presence of a 60 • dislocation at the Ge/Si interface lowers the barrier for nucleation of a second (the "complementary one") 60 • dislocation with Burgers vector such that the superimposition of the two generates a 90 • dislocation. In order for this reaction to take place through simple dislocation gliding, however, the original dislocation must lay on the glide plane of the complementary one. While this can happen in some cases, the most probable configuration leads to close pairs of 60 • dislocations, clearly observed in low-temperature experiments 27,32 . How such pairs evolve into a 90 • dislocation upon annealing is not yet clarified. Similarly, no clear answer to question (ii)-regarding annealing-induced ordering-is present in the literature. The possible role of vacancies in influencing gliding/climbing and lateral ordering was invoked 21 but only at the speculative level. As reaching defect uniformity is an important requirement for applications, we find it of utmost importance to provide a solid atomic-scale explanation.
In this work we shall clarify both the mechanism of 90 • dislocation formation from pairs of complementary 60 • dislocations and of lateral ordering at the Ge/Si interface by suitable molecular dynamics (MD) simulations, highlighting the key role played by vacancies in influencing defects motion and the importance of the growth/ post-growth conditions.

Methods
All simulation results presented in this work were obtained in the framework of classical molecular dynamics as implemented in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 36 . The Tersoff potential 37 was used to describe interatomic interactions. While obviously not exact, such an approach has been widely exploited to model defects in Ge/Si systems 32,[38][39][40] , yielding results compatible with experimental evidence. Simulations presented in the present article were run in the canonical ensemble, using a Nose-Hoover thermostat 41 and setting the time step to 1 fs. The latter was calibrated in a set of dedicated microcanonical simulations where relative fluctuations in total energy lower than 10 −5 were observed even at the highest ( T = 1400 K ) sampled temperature.
Orthogonal simulation cells, oriented along the directions determined by the vectors x = [110] , y = [110] , z = [001] were used, and periodic boundary conditions (PBCs) were applied in the x and y directions. To model a strained epitaxial layer, we created a simulation cell composed of 32 monolayers of Ge atoms on top of other 32 monolayers of silicon atoms along the z = [001] direction. The lattice parameter was set as the one of pure silicon. In this way germanium atoms were strained in the x, y directions ε xx = ε yy ≈ 4% . To simulate the effect of a thick bulk substrate below the Ge layer, the bottom three layers of silicon atoms in the simulation cell were kept fixed to bulk positions. The top Ge surface was reconstructed 2 × 1 upon forming dimers accordingly. We performed simulations of dislocation arrays by using three different cell dimensions ( 9.987/14.980/19.973 × 1.536 × 8.530 nm 3 ) depending on the explored strain conditions. The cells are made up by replicating the perfect diamond lattice unit cell 26/39/52 × 4 × 16 times.
Atomic trajectories were analyzed using the Open Visualization Tool (OVITO) software 42 , which allows for a very convenient identification of defects in crystals, e.g., providing dislocation lines along with the associated Burgers vectors.
Dislocation modeling. For the sake of simplicity, dislocations were introduced in the system as straight segments running along the y direction. Starting from a more realistic (half)loop configuration would have required a larger cell in the y direction, making it impossible to reach the remarkable time scales simulated when describing climb (see the "Results" section). Instead, we could set the cell length along y as relatively small ( 1.536 nm ). Shortening the dislocation line further would be critical as it would preclude the formation of kinks and jogs, allowing for the natural motion of dislocations. On the other hand, despite being minimal, we verified that the simulation cell was indeed large enough along y: we repeated a subset of the simulations using a three-time larger cell along y without observing any significant dependence of the results on the dimensions (see Supplementary Information online).
Dislocation lines were inserted in the simulation cell by shifting each atom by the displacement field vectors calculated by dislocation modeling in the framework of the linear elasticity theory 43 . The obtained configurations were geometrically optimized using a conjugated gradient minimization algorithm. We recall that, for an arbitrary dislocation segment, the field components can be easily calculated by rotating the reference system and decomposing the dislocation Burgers vector into a screw ( b ) and an edge component ( b ⊥ ). At this, ẑ axis has to be oriented along the dislocation line ξ , x axis perpendicular to ξ in the glide plane, and ŷ axis oriented following the right hand rule, respectively. In this case, the components of the displacement field vectors are written as follows:  (1) gives the displacement field associated with an isolated dislocation in bulk. Before applying it to the atoms in the simulation cell we found it convenient to also add the contribution of the image dislocation 43 , introduced by the presence of the free surface. We also added the additional contribution induced by the presence of PBCs along the x direction 44 creating a dislocation array out of a single defect.
Simulation of vacancy-induced dislocation motion. Dislocation climbing requires vacancies to migrate to the dislocation core along the whole dislocation line 43 . Brute-force runs, where both typical experimental temperatures and physically-sound vacancy concentrations are considered, do not show any significant evolution due to the timescale limitations imposed by MD. Even if one disregards experimental conditions, observing climbing along the full dislocation line with MD is not trivial. One can attempt to speed up the motion by increasing the simulation temperature, raising the vacancies' diffusion coefficient. However, there exists a limit beyond which excessive disorder is introduced in the system (we used up to T = 1400 K , high but well below the Tersoff-potential Ge melting temperature 45 ). Instead of manipulating the diffusion coefficient, other approaches can be used. For instance, Li et al. pre-positioned an entire row of vacancies on the dislocation path, therefore forcing climbing motion 46 once the dislocation reached them. A step forward would require to simulate also the motion of vacancies towards the dislocations. To achieve this within reasonable computational time, it is required to reduce the path a vacancy must travel to reach the dislocation core. A possible way to accomplish this is to consider very high vacancy densities. If vacancies are placed too close to each other, however, the formation of slowly-diffusing vacancy complexes precludes observation of migration to the dislocation core. After attempting several strategies along the lines traced here-above, we came up with a procedure that we find particularly appealing since the system is left sufficiently free to evolve, and climbing of the entire dislocation line can be observed. The procedure is the following one: we start by introducing a dislocation (an extension to two dislocations is also exploited in the "Results" section) in the simulation cell. We then remove one atom of the crystal randomly. We start the MD simulation and track the vacancy motion on the fly by performing a Wigner Seitz Defect Analysis using the OVITO software. If the vacancy ends up at the dislocation core or at the top free surface, we then randomly insert a second vacancy and start tracking its motion. Otherwise, we keep on simulating the system dynamics. The procedure is iterated until the desired number of vacancies are inserted. As demonstrated in the "Results" section, the procedure works and allows to observe the vacancy-induced motion of a full dislocation line. In principle, we could have randomly created vacancies within the whole simulation cell. However, to decrease the time needed for the vacancy to reach the dislocation core, we defined a smaller "vacancy-spawn" regionwhere the point defect was created-composed of the atoms distant less than 3 nm from the dislocation core. The 3 nm size of the "vacancy-spawn" region is relevant only for speeding up the simulations. Considering a "spawn-region" large as the whole cell along x did not change the ratio of vacancies that effectively migrate to the dislocation core. In addition, such a region was restricted to be at a distance from the upper free surface larger than the one of the dislocation core as vacancies created between the dislocation core and the free surface were observed to annihilate at the latter irreversibly.

Results and discussion
Gliding of pairs of 60 • misfit dislocation. As already stressed in the introduction, induced nucleation of a second complementary 60 • dislocation by a first 60 • dislocation at the Ge/Si interface has been recognized as the primary cause for the creation of 90 • dislocation and of 60 • dislocation pairs 32,47 . As no atomic-scale description of further evolution under annealing is present in the literature, we shall start our analysis precisely from these pairs. Following Ref. 32 , the typical initial condition leading to their formation is given by having one 60 • dislocation-called "A" in Fig. 1-at the Ge/Si interface and a second one-"B" in Fig. 1-nucleated close to the free surface. If the position of the core of dislocation "A" belongs to the glide plane of dislocation "B" (this glide plane is often called mirror glide plane, and is indicated by a red dashed line in Fig. 1) then simple gliding of "B" towards "A" can lead to the formation of a 90 • dislocation. For this case, annealing is not needed to promote the existence of such a defect. A close pair of 60 • dislocations is formed when, instead, dislocation "B" glide plane intersects dislocation "A" glide plane not at the interface. Three representative initial configurations are displayed in Fig. 1, together with the evolution predicted by molecular dynamics simulations at 1000 K . Dislocation "A", with Burgers vector � b = a/2[011] is placed at the interface (the dislocation core is highlighted in green), dislocation "B", with Burgers vector � b = a/2[101] , is placed 1.5 nm below the free surface (the dislocation core is highlighted in red). The glide planes of the two dislocations are misaligned (they do not intersect at the interface); therefore, they cannot form a 90 • dislocation at the interface. The distance d ML , between the glide plane of the second dislocation (at the interface) and the mirror-like glide plane, is used as a parameter. Three cell sizes along the x direction were used to simulate different relaxation of the Ge film: ≈ 10 nm-corresponding to full relaxation of the ∼ 4% Ge/Si lattice mismatch in the presence of a 90 • edge dislocation at the interface per cell, ≈ 15 nm-i.e., 75% relaxation with a 90 • edge dislocation at the interface, ≈ 20 nm-i.e., half relaxation with a 90 • edge dislocation at the interface. www.nature.com/scientificreports/ A wide range of relative arrangements of two complementary dislocations are possible and was observed in low-T experiments. Still, the results of their evolution can be grouped in a few categories only: an edge dislocation at the Si/Ge interface (see Brunetto et al. 32 ), an edge dislocation few nanometers above the Si/Ge interface, a 60 • dislocation "floating" few nanometers above the Si/Ge interface in the vicinity of its complementary counterpart, and two complementary 60 • dislocations both at the interface. MD simulations were able to fully capture all these behavior, as seen in Fig To explain the results of the here-above reported molecular dynamics simulations, we find it helpful to show, in Fig. 2, the force field experienced by dislocation "B" placed in different positions while keeping dislocation "A" at the Ge/Si interface. Such force field can be conveniently computed based on elasticity theory 48 , including periodic boundary conditions and the effect of the top free surface. Figure 2 allows for an immediate interpretation of the various configurations obtained at the end of the molecular dynamics simulations of Fig. 1. Note that, differently from molecular dynamics, the analytical calculation assumes 48 an isotropic medium. However, it is evident from the comparison between Fig. 2 and the final configurations of Fig. 1 that the qualitative behavior found by MD and predicted by the analytical model corresponds.
The force field displayed in Fig. 2a is computed for an array of dislocations with a periodicity of 10 nm (ideal periodicity of 90 • dislocations to relax pure Ge on Si). The one shown in Fig. 2b, instead, corresponds to a 20 nm periodicity. The force vectors are colored in red/blue to show the intensity and the orientation of the glide component of the force in a specific position. A saturated red means a strong glide force towards the Si/Ge interface, saturated blue means a strong glide force towards the free surface, pale/white hues indicate a small/zero force. In Fig. 2 a star indicates where dislocation "B" will stop its glide for some selected glide planes. The formation of an edge dislocation at the Si/Ge interface is achieved when dislocation "B" nucleates and glide along a plane which www.nature.com/scientificreports/ is exactly the mirror-like glide plane (double solid line in Fig. 2). In this case, the force acting on the dislocation is almost aligned to the glide plane resulting in an easy glide of the dislocation. An edge dislocation above the Si/Ge interface can be formed when dislocation "B" nucleates slightly apart from the mirror-like glide plane. Consider for reference the single solid line in Fig. 2; in this case, dislocation "B" can move on a glide plane which is 1.5 nm away from the mirror-like glide plane. The projection of the force acting on dislocation "B" on its glide plane (i.e., the glide component of the force) changes in sign a few nanometers above the Si/Ge interface (see red/blue color-changing). Therefore the dislocation, in this case, will glide towards the Si but will stop in the vicinity of its complementary dislocation (a star indicates the position where the dislocation will stop in Fig. 2). It clearly corresponds to what happen in the molecular dynamics simulation of Fig. 1a-d. At that height dislocation "B" can attract dislocation "A", making it gliding upward to merge and form a 90 • dislocation at their junction (Fig.1d).
A qualitatively different result is obtained by placing dislocation "B" further away from the mirror-like glide plane. Indeed, when the distance between the glide plane of "B" and the mirror-like glide plane ( d ML ) is 3.5 nm (dashed black line in the Fig. 2a) we can notice that at a certain height the force acting on the dislocation is directed for the most outside its glide plane. This can result in an evolution which is analogous to that reproduced via molecular dynamics and shown in Fig. 1e-h: dislocation "B" migrate from the free surface towards the Si/Ge interface, but in the middle of the film it stops in a local minimum position because of the small glide component of the Peach and Koehler force. Dislocation "B" remains floating in the Ge film, too far away from dislocation "A" to make it glide towards it.
This behavior can be expected to change in the case of a much higher strain. Indeed using a periodicity of 20 nm instead of 10 nm , as shown in Fig. 2b, the glide component of the Peach and Koehler force is much higher, and dislocation "B" is pushed towards the Si/Ge interface with almost the same force intensity throughout all the Ge film. In this latter case dislocation "B" will reach the interface, exactly as reproduced via the molecular dynamics simulations illustrated in Fig. 1i-l. We have checked that in the intermediate strain condition ( 15 nm periodicity), the final dislocation configuration is similar to those reported in the 10 nm map.
In many low-temperature growth experiments, proofs of the presence of these particular arrangements of complementary dislocations are commonly observed 27,32 , but when the successive annealing at high temperature is performed, most of the defects are found at the interface, they are ordered, and most of them are 90 • edge Climbing of 90 • misfit dislocation. We analyzed the role of the climbing mechanism in the formation and ordering of 90 • edge dislocations starting from the configuration found as the final state reached by glide motion. In particular, we first considered the configuration of dislocations in Fig. 1d-i.e. a 90 • Lomer dislocation formed above the Si/Ge interface. Climbing of an edge dislocation is the simplest case of dislocation climbing motion and is often used in the literature as a prototypical example 43 . We have already described in the "Methods" section how we triggered climbing by inserting vacancies one at the time, removing atoms of the crystal randomly within a "vacancy spawn region, " before simulating the system evolution at T = 1400 K and checking whether the conditions to insert a new vacancy are met. Here we shall illustrate the procedure further by analyzing the actual simulation where climbing of the 90 • edge dislocation is tackled (Fig. 3). As shown in Fig. 3e, the first vacancy, which is represented by the large blue sphere, is observed to diffuse through the Ge film (blue trajectory in the same Figure) while the edge dislocation stays immobile. After few hundred ps, the vacancy reaches the dislocation core. At this point another vacancy is inserted (red sphere in Fig. 3e) and the system is evolved until also the second vacancy reaches the dislocation core (see Fig. 3f, this time its trajectory is highlighted in red). The simulation is iterated, and a new vacancy is created in the "vacancy spawn" region every time there are no vacancies in that region, i.e., the vacancy has reached the dislocation core or is lost towards the free surface. The third vacancy (green sphere) reaches in turn the dislocation core (Fig. 3g), and thereafter, in few hundreds of ps, is followed by the fourth vacancy (yellow sphere, Fig. 3h). With four vacancies, there is a section of the dislocation line that is one layer below its original position, i.e., a jog has formed. A fifth vacancy is inserted and it migrate to the dislocation core (tile sphere, Fig. 3i). The successive vacancy is the sixth one, and it also reaches the dislocation core that now is almost completely one glide plane below its starting position (purple sphere, Fig. 3j). The seventh vacancy (gray sphere, Fig. 3k) does the same, and finally the last vacancy migrates to the dislocation core, achieving the complete climb of the dislocation (orange sphere, Fig. 3l). The starting and final configuration are shown them Fig. 3a, c, respectively. A black arrow show the displacement of the dislocation core from its starting position. The top views of the dislocation core with and without vacancies are shown in Fig. 3d, b, respectively. The analysis of vacancy trajectories clearly shows that the motion of a single vacancy is random-like at the beginning, while in the proximity of the dislocation core, the attraction between the point and the linear defect leads to the irreversible capture of the former at the core of the latter. A video animation of the climb of an edge dislocation on a three times larger cell can be found as Supplementary Video S1 and Supplementary Fig. S1  www.nature.com/scientificreports/ online. It is also evident that vacancies are attracted by the compressive deformation below the dislocation core. Indeed, a relevant trend, which can be easily observed by considering trajectory lines, is that once vacancies are close to the dislocation core, they start to move towards the dislocation preferentially along {111} planes. Note that the dislocation's climb movement turns to be a step towards the interface along a 111 direction.
Climbing motion involving pairs of 60 • misfit dislocation. The configurations that we have found by pure gliding (Fig. 1) and considered for the formation of purely 90 • edge dislocation at the interface were grouped into three classes. The first set was made up by an edge dislocation a few nanometers above the interface, and its evolution was described in the preceding section. Then two other sets were possible, both composed of complementary 60 • misfit dislocation pairs. Being clarified that, when we follow the evolution from the first configuration, 90 • dislocation can be found at the Si/Ge interface because of climbing motion activated by vacancies; it remains to be understood if the same stems for pairs of complementary 60 • misfit dislocations. Two configurations are possible: two 60 • complementary dislocations spatially separated at the interface (see Fig. 1i-l for reference) or a 60 • dislocation blocked in the film above the interface in the vicinity of its complementary pair (such as in Fig. 1e-h). After having demonstrated that the introduction of one vacancy at a time in the simulation cell allows one to conveniently simulate climbing of an edge dislocation, we extended the same procedure for investigating vacancy-induced motion of 60 • pairs. In particular, we started by investigating the evolution of a 60 • dislocation blocked a few nanometers above its complementary pair at the interface. Dislocations are positioned as shown in Fig. 1h, and we will refer to dislocation "A" and "B" in the same way. MD simulations (not reported here) showed us that dislocation "B" remains immobile (stuck in a local minimum) until enough vacancies make it climb by a single glide plane. After that, the dislocation glides again, eventually reaching the interface. Notably, the response of a 60 • misfit dislocation, which is glissile, to the interaction with vacancies is not a pure climb motion. The formation of a 90 • Lomer at the interface can still be achieved by the motion of two complementary 60 • dislocations alongside the Si/Ge interface. The situation where two complementary 60 • misfit dislocations are located nearby at the interface has been revealed in many experiments, in particular before annealing 23,47 . Annealing can cause vacancy migration towards dislocations, making them climb one towards the other, eventually leading to dislocation joining, forming an edge dislocation. To check this hypothesis, we studied the interaction of vacancies with two complementary 60 • dislocations already at the Si/Ge interface as in Fig. 1k, separated by a distance of 3 nm . We inserted a vacancy at a time, randomly deleting an atom in a region below the dislocation cores, excluding the cores themselves by the region (excluding a cylindrical area with axis along the y direction, a radius equal to 0.7 nm , centered in the middle point of the two dislocations).
In Fig. 4 we show snapshots of MD simulations of the evolution of vacancies in the presence of the two complementary 60 • dislocations. In particular, in Fig. 4a the initial position of the two dislocations is shown. The simulated time and the total number of vacancies inserted so far in the cell are shown in the insets. Each dislocation core is highlighted in a different color, we will refer to them following the same convention used in Fig. 1, dislocation "A" for the one with Burgers vector � b = a/2[011] (colored in green) and dislocation "B" for the one with Burgers vector � b = a/2[101] (colored in red). Dashed lines, colored accordingly, also show their glide planes. After 13 vacancies have been inserted in the cell and a simulated time of more than 14 ns , dislocation "A" climbs, moving outside its original glide plane (the dotted green line) towards the next glide plane (dashed green line), see Fig. 4b. Now dislocation "A" has climbed inside the Si/Ge interface, and therefore it glides along its glide planes back to the Ge film, moving towards dislocation "B" (Fig. 4c). After that, in the successive two ns, dislocation "A" climbs and glide again, as displayed by the arrows in Fig. 4d. The former glide planes are shown for reference as dotted green lines, with the thinner being the original one. Also dislocation "B" climbs, as shown in Fig. 4e, reducing the distance between the two dislocations. Actually, at such limited distance, the attraction between the two dislocations becomes so intense that dislocation "B" glides along its glide plane and reaches the same height above the interface of dislocation "A. " Still, some more vacancies are needed to make the two dislocations climb and react. Enough vacancies eventually arrive at the dislocation cores after 41.2 ns , as shown in Fig. 4g (old glide planes are shown by dotted lines for reference). The glide plane of the newly formed 90 • edge dislocation is 2 mono-layers above the Si/Ge interface (see Fig. 4h). Additional migration of the newly formed 90 • edge dislocation (which core is highlighted in yellow) happens. It climbs to the Si/Ge interface adding more vacancies in the cell, as shown in the last snapshot (Fig. 4i). The new glide plane is shown by a dashed yellow line, and the former one by a dotted yellow line, an arrow indicates the displacement of the edge dislocation.
In Fig. 4 we have chosen to present a section of the whole-cell taken by slicing the cell perpendicularly to the dislocation lines to ease the view. In effect, there is some disorder at every step of the present simulation because we are simulating the evolution of 60 • dislocations that are not sessile, as in the case of an edge dislocation. Therefore dislocation segments move easily and not as a whole. As a result, we observed kink formation along the dislocation lines. In order to obtain as final configuration a straight perfect edge dislocation at the interface long as the whole simulation cell, we had to add a 50ns equilibration phase during which the frequency of vacancy insertion was decreased. While not the main aim of this work, we notice how the present simulations, see Fig. 4i, naturally predict the presence of Ge/Si intermixing around dislocation cores and in the surrounding areas, as observed in experiments and justified theoretically based on elastic-energy minimization 23,49 . In such references, the actual microscopic kinetic mechanism was not described. The present simulations seem to show that mixing can be naturally produced by the motion of the vacancies towards dislocation cores. Our simulations also partially reproduced this mixing of Ge atoms in the substrate and vice-versa (note Ge atoms colored in blue in the Si region and vice-versa in Fig. 4f-i).
Glide of 90 • misfit dislocation. In the former two sections, we investigated the role of vacancies in dislocation climb. As expected, simulations showed that vacancies are attracted by dislocations, eventually pro- www.nature.com/scientificreports/ ducing their climbing motion. As a result, edge dislocations from the Ge film can migrate to the Si/Ge interface, and pairs of complementary 60 • dislocations can migrate alongside the interface and there join to form edge dislocations. Based on our results, we can therefore provide the following explanation of the experimental findings 21,29,32 : during low-T growth, a large density of vacancies is introduced in the system as a direct consequence of out-of-equilibrium growth. However, their mobility is limited once again due to the low T. The situation changes when a second stage involving either annealing or high-T growth follows. Vacancies can now quickly diffuse, enabling the interaction of complementary 60 • misfit dislocation, originally grouped into pairs, eventually forming edge dislocations.
Having provided an answer to the first open question regarding the annealing-induced abundance of 90 • versus 60 • dislocations, we can now address the second, i.e., how such edge dislocations can move laterally, producing an ordered distribution as observed experimentally 21,29 .
To answer, we considered the situation where an array of non-equidistant 90 • is present at the Si/Ge interface. We considered a simulation cell 20 nm wide. We inserted two edge dislocations distanced by only 3 nm (as we use PBCs, this actually corresponds to an array of dislocations where each defect has a distance of 3 nm with the closest defect and of 17 nm with the next one) checking whether mechanism leading to an increase of such distance towards the 10 nm ideal value was occurring.
A Lomer dislocation has Burgers vector � b = a/2[110] or � b = a/2[110] , and its glide plane is the (001) plane, that coincides with the Ge/Si interface. That plane does not belong to the primary slip system of the zinc-blende/ diamond structure. Hence, the glide motion of such a dislocation is not as easy as for 60 • dislocations. Even if vacancies are not necessarily required for the glide motion of edge dislocations, we will show that their presence is of fundamental importance in accelerating dislocation movement, enabling Lomer dislocation glide on a reduced timescale.
Indeed, we ran an MD simulation at 1400 K without any vacancy in the Ge film, and we observed a single gliding event (i.e., lateral motion leading to an increase in lateral distance) after an extremely long simulation time ( ∼ 90 ns ). At a lower T of 1200 K, no events at all were witnessed in as long as 100 ns , the maximum time we reached with our simulations. We conclude that lateral motion can occur but on quite long time scales. Interestingly, however, we noticed that a single vacancy could lead to a noticeable increase in gliding velocity. In Fig. 5, indeed, we show the results of a simulation where a single vacancy is inserted a few lattice sites to the right of Lomer dislocation "C", in Fig. 5a. The vacancy, represented by a purple sphere, is attracted by the dislocation and reaches its core in a relatively short timescale (Fig. 5b). No vacancy is present near the other edge dislocation www.nature.com/scientificreports/ "D. " The presence of a vacancy in the dislocation "C" core enables the formation of a kink as shown in Fig. 5c (the top view of the defect, together with the corresponding dislocation line, is provided).
Once the kink has formed, the completion of the glide step is relatively fast; indeed, only ≈ 100 ps are needed to move the whole dislocation line (see Fig. 5d-e). The vacancy is still around the dislocation core so that it can enhance the glide motion of dislocation "C" again. Indeed, after ≈ 50 ns dislocation "C" glides on the right again as shown in Fig. 5f (the black arrows are inserted for reference). A third and a fourth glide jumps to the right are observed in the next ≈ 30 ns (Fig. 5g-h). Notably, we also observe dislocation "D" glide on the left, but after 92.3 ns of simulation time. The video animation of the activation of the lateral glide of a 90 • edge dislocation via vacancy absorption on a three times larger cell along y can be found as Supplementary Video S2 and Supplementary Fig. S2 online.
In the presence of the vacancy, a single glide step was also observed in a simulation run at 1200 K.

Conclusions
In this work, we exploited molecular dynamics simulations to shed light on the atomic-scale mechanisms leading to the experimentally observed formation of a laterally-ordered array of 90 • dislocations at the interface of Ge/ Si(001) films [21][22][23]29 . Notably, such arrays are typically obtained in experiments where a first low-T deposition stage is followed by a rise in temperature 21,29 . In the former stage, pairs of 60 • dislocations are formed by pure gliding motion, which can occur even in the absence of fast-moving vacancies. And, indeed, we could observe the formation of such pairs in MD simulations where no point defects at all were present. Explaining the complete transition of the pairs into 90 • dislocations required climbing motion to be activated at the higher temperature following the low-T stage. After having introduced a convenient procedure to introduce vacancies in the simulation cell, we were able to observe full vacancy-induced dislocation climbing leading to the transformation of the 60 • pairs into 90 • dislocations and/or climbing of 90 • dislocations at the Ge/Si(001) interface. At this point, we tackled lateral ordering of 90 • dislocations at the interface. While lateral gliding motion was also observed in the absence of point defects, our simulations revealed how even the presence of a single vacancy could significantly speed up the process, ultimately leading to an ordered array. With this respect, we find it interesting to point out that the final stage of lateral ordering can also be reached without imaging an excess of vacancies in the system, at variance with the transformation of 60 • into 90 • in the likely situation where nucleation of the second 60 • dislocation does not occur on the mirror glide plane.
In conclusion, our simulations reproduced experimental findings and, in doing so, helped in understanding the different phases needed to achieve the limit of a unimodal distribution of laterally-ordered dislocations in www.nature.com/scientificreports/ Ge/Si(001), highlighting when the role played by vacancies is crucial and when it is not. This can help to devise new growth strategies where alternating close-to-equilibrium and far-from-equilibrium conditions are exploited to influence the dislocation distribution.